A dynamic calcium-force relationship model for sag behavior in fast skeletal muscle

In vitro studies using isolated or skinned muscle fibers suggest that the sigmoidal relationship between the intracellular calcium concentration and force production may depend upon muscle type and activity. The goal of this study was to investigate whether and how the calcium-force relationship changes during force production under physiological conditions of muscle excitation and length in fast skeletal muscles. A computational framework was developed to identify the dynamic variation in the calcium-force relationship during force generation over a full physiological range of stimulation frequencies and muscle lengths in cat gastrocnemius muscles. In contrast to the situation in slow muscles such as the soleus, the calcium concentration for the half-maximal force needed to drift rightward to reproduce the progressive force decline, or sag behavior, observed during unfused isometric contractions at the intermediate length under low-frequency stimulation (i.e., 20 Hz). The slope at the calcium concentration for the half-maximal force was required to drift upward for force enhancement during unfused isometric contractions at the intermediate length under high-frequency stimulation (i.e., 40 Hz). The slope variation in the calcium–force relationship played a crucial role in shaping sag behavior across different muscle lengths. The muscle model with dynamic variations in the calcium-force relationship also accounted for the length-force and velocity-force properties measured under full excitation. These results imply that the calcium sensitivity and cooperativity of force-inducing crossbridge formation between actin and myosin filaments may be operationally altered in accordance with the mode of neural excitation and muscle movement in intact fast muscles.


Introduction
Skeletal muscles are specialized tissues that can contract to generate force in response to neural signals from spinal motoneurons. It is commonly thought that contractile forces are induced from crossbridge formation between actin and myosin filaments in the sarcomere [1][2][3][4], which is activated by calcium release from the sarcoplasmic reticulum according to neural excitation [5]. Thus, the transfer function of force inducing crossbridge formation across different excitation levels has been characterized by measuring steady-state force responses at various constant calcium concentration levels [6].
Isolated or skinned muscle fiber experiments on the calcium sensitivity of muscle force have shown larger force production at higher calcium concentrations compared to the lower constant calcium concentration levels at steady state [7]. This sigmoid steady-state calciumforce relationship implies the cooperativity and saturation of crossbridge formation as the level of muscle excitation increases [8][9][10]. The steady-state calcium-force relationship has also been shown to be specific to the type of muscle fiber (slow-twitch versus fast-twitch in [11]) in mammals [12][13][14][15][16] and humans [17]. Recent in vitro studies have further demonstrated that the shape of the calcium-force relationship can be considerably altered by inorganic phosphate molecules known to accumulate during force production [18] or variation in muscle length [7], presumably owing to changes in actin-myosin interaction [19]. However, it remains elusive whether and how the calcium-force relationship varies in intact muscle fibers under physiological conditions of muscle excitation and length.
Computational studies have suggested that the calcium-force relationship might differ depending on the type of muscle fibers (slow-twitch versus fast-twitch) under physiological conditions. Slow-twitch muscles have been shown to produce a force whose peak is maintained during isometric contractions at any level of stimulation frequency [20]. These isometric forces could be reproduced for a full physiological range of stimulation frequencies with no changes in the calcium-force relationship in computational models [21]. In contrast, fasttwitch muscles have been shown to produce a force that declines rapidly within half a second after its peak during isometric contractions under unfused stimulation conditions [22]. The impaired crossbridge formation by the inhibitory effects of inorganic phosphate accumulation has been computationally and experimentally demonstrated as a potential mechanism for this rapid force decline (or sag) during continuous submaximal contractions [23]. On the other hand, the decline in free calcium ion concentration in the sarcoplasm by the reduction in calcium release from the sarcoplasmic reticulum [24] and the accelerated uptake of calcium into the sarcoplasmic reticulum [25] have been proposed to be involved in the final phase of muscle fatigue, showing accelerated force reduction several seconds after the onset of continuous maximal contraction [24,26].
Here, we focus on the dynamic variation in the calcium-force relationship during the rapid force decline (or sag) in fast skeletal muscles before the final phase of muscle fatigue. It is hypothesized that both the midpoint and the slope of the sigmoid calcium-force relationship dynamically change for force production under physiological conditions. Because of current experimental limitations, a computational approach is developed to test our hypothesis. The comparison between experimental and simulation results shows that dynamic variation in the calcium-force relationship is required to replicate the force responses of the cat medial gastrocnemius over the full physiological range of stimulation frequencies and muscle lengths. This study provides insights into the dynamic properties of calcium-sensitive crossbridge formation under physiological conditions and a computational framework for realistic modeling of complex muscle behaviors.

Ethics statement
All procedures for collecting data from cat medial gastrocnemius (MG) were approved by the Animal Care Committee at Northwestern University.

In situ experiments
The properties of average active force production by single sarcomeres were characterized by measuring the active force of the whole muscle in decerebrated cats [3]. Details of the surgical preparations and experimental protocols have been previously reported [27]. Briefly, adult cats (CAT14 for model development and CAT12 for prediction confirmation) were anesthetized, the medial gastrocnemius was exposed, and the distal tendon was attached to a computer-controlled puller. All distal hindlimb nerves except the medial gastrocnemius nerve were disconnected, only ipsilateral dorsal roots from the fourth lumbar spinal nerve (L4) to the second sacral spinal nerve (S2) were transected to eliminate sensory feedback from the medial gastrocnemius, and decerebration was performed at the midbrain. Hindlimb and core temperatures were maintained within physiological limits. Electrical impulses with a width of 0.1 ms were applied to control muscle force either using stainless steel wires in the proximal and distal portions of the muscle belly or via hook electrodes on the ventral roots. A supra axial current (50-100% above that required to elicit full recruitment) was used to produce repeatable and consistent forces during all trials. The length of the muscle-tendon unit was controlled by a position servo system and measured by a linear variable differential transformer attached to the puller shaft. Muscle forces were measured by a strain gauge-based transducer in series with the shaft. Raw data for muscle force and length were smoothed using the function ("smooth" with options of "0.01" for span and "loess" for method) built in MATLAB software (version R2012a, MathWorks, Natick, MA). To isolate the active force response, the force data measured without the stimulation (i.e., passive response) were subtracted from those measured with the stimulation (i.e., total response) under the same condition of muscle-tendon length. The passive and active trials were separated by a 30 s rest period, and all active trials were also separated by at least one minute to minimize fatigue.
intermediate muscle-tendon length (X m,0.5 , 5 mm for CAT14 and 3.64 mm for CAT12) between the physiologically minimal (X m,0 , 0 mm for CAT14 and -0.55 mm for CAT12) and maximal (X m,1 , 10 mm for CAT14 and 7.9 mm for CAT12) lengths that have been observed during locomotor-like movements in cats [28][29][30]. The representative types of sag behavior were reproduced under constant stimulation frequencies of 20 and 30 Hz at X m,0.5 for both CAT14 and CAT12. The length and velocity dependence of the dynamics of the calcium-force relationship was investigated during shortening and lengthening at constant velocities and step lengthening over a very short period of time under full excitation (100 Hz) for both CAT14 and CAT12.

Modular modeling framework
In this study, we used a modular modeling approach that allowed us to divide the complex process of force generation into subprocesses that can be separately modeled and validated by relevant experimental data. A modular model of the muscle-tendon complex experimentally validated for the adult cat soleus (predominantly slow twitch fibers [20]) in in situ conditions was adapted to replicate the force production of the adult cat medial gastrocnemius (predominantly fast twitch fibers [31]). The derivation and description of model equations for the model of slow skeletal muscles has been fully presented in a previous study [21]. We employed Hill-type mechanics comprising a contractile and an elastic element in series [32] because it is simple and effective for modular modeling of the essential processes that intact muscle undergoes for contractile force production under physiological conditions (Fig 1). The first module (i.e., Module 1) transforms the electrical signals coming from either spinal motoneurons or electrical stimulation into the concentration of free calcium (Ca) within the sarcoplasm (SP). The concentration of sarcoplasmic calcium was determined by multiple factors, including the reaction of calcium (Ca SR ) and calsequestrin (CS) in the sarcoplasmic reticulum (SR), calcium release (R) and uptake (U) through the membrane of the SR, calcium buffering by free proteins (B), and calcium-troponin binding (CaT) modulated by the muscle-tendon length (X m ) and muscle activation level (Ã) in terms of the constant CaT. The system equations for Module1 are as follows: where CS 0 , B 0 and T 0 indicate the concentration of total calsequestrin in SR, total free calciumbuffering proteins, and total troponin in SP, respectively, and K1-K6 are the rate constants for chemical reactions between Ca SR , Ca SR CS, Ca, B, T, CaB and CaT.
The action of the calcium pump in the SR was mathematically modeled as R and U in the absence of an inhibitor for pump action [33], where R max , τ 1 , τ 2 , U max , and K indicate the maximum calcium release rate from the SR, the time constant of calcium release increase, the time constant of calcium release decrease, the maximal calcium uptake rate, and the association constant for calcium binding to pump, respectively.
To reflect the length-dependent changes in the calcium-force relationship [34], the default value for the forward rate constant (i.e., K5 i ) of the calcium-troponin reaction was made as a and muscle mechanics (F) induced by contractile elements (CEs) and serial elastic elements (SEs). X CE and X m show the length of the contractile element and entire muscle-tendon unit. C1 and 1/C2 indicate the constant concentration of calcium-binding troponin relative to the total troponin concentration (CaT/T 0 ) for the half-maximal level of steady-state muscle activation (Ã~1) and the slope of the CaT/ T 0 -Ã~1 curve at C1. B. Relationship of steady-state C1 (C1 1 ) and C2 (C2 1 ) with the concentration of CaT/T 0 (upper panel) and time constants (τ C1 and τ C2 ) for the C1 and C2 dynamics (bottom panel). C1n1 and C2n1 limit the saturation level, and C1n2 and C2n2 indicate a constant concentration of CaT/T 0 for the midpoint between the minimum and maximum in the C1 1 and C2 1 curves. function of X m to match twitch responses at different X m values, as suggested in a previous study [35], In addition, the affinity increase of troponin for calcium by crossbridge attachments [36] was reflected by varying the default value for the backward rate constant (i.e., K6 i ) inversely tõ A in a similar manner as reported in a previous study [33], The second module (i.e., Module 2) transforms the concentration of CaT relative to the total troponin concentration (T 0 ) into the level of muscle activation (A), representing a fraction of actomyosin crossbriges available for force generation. A was determined through exponentiation ofÃ with a single exponent (α) to reflect the degradation of muscle activation for the fluctuation in CaT under physiological conditions. The system equations for Module 2 are as follows: The dynamics ofÃ were mathematically represented, The third module (i.e., Module 3) transforms the muscle activation level (A) resulting from crossbridge formation into the force (F) based on Hill-type mechanics. F was determined via multiplication of the stiffness (K SE ) of serial elastic elements (SE), including extracellular (i.e., tendons and aponeuroses) and intracellular (i.e., titin and myofilaments) components, with the variation in SE length reflecting the length-and velocity-tension properties measured under full excitation. The SE length variation was calculated by subtracting a change in the contractile element (ΔX CE ) from a change in the whole muscle-tendon unit (ΔX m ) length. The length-tension relationship was proportionally adjusted to the level of muscle activation (A) resulting from neural excitation or electrical stimulation [37]. The tension axis of the velocitytension curve was scaled by the length-tension relationship to reflect the overlap of sliding filaments determining the number of force-generating crossbridges [38]. The system equations for Module 3 are as follows: where P 0.5 is the peak force at the intermediate length between the physiologically minimal and maximal length under full excitation in the isometric condition and K SE is the stiffness of the serial element normalized by P 0.5 .

Dynamics of the calcium-force relationship
In module 2 of our muscle-tendon model, the calcium-force relationship was phenomenologically represented due to the lack of experimental data for underlying molecular mechanisms, and the purpose of this study was to investigate the dynamics of calcium sensitivity and cooperativity for force production under physiological conditions. The calcium-force relationship was modulated by dynamically varying the two factors that shape the sigmoid steady-state calcium-activation (CaT/T 0 -Ã 1 ) relationship (see Fig 1 and Eq (6) for the graphical and mathematical description). One factor (C1) representing calcium sensitivity was the concentration of calcium-binding troponin to reach half of the maximal muscle activation. The other (1/C2) representing cooperativity was for the slope of the calcium-activation relationship at C1. C1 and C2 were mathematically formulated as a sigmoidal function that considers the concentration of calcium-binding troponin relative to the total troponin concentration (i.e., CaT/T 0 ). The equations for the dynamics of C1 and C2 are as follows: where X indicates 1 for C1 and 2 for C2; CXn1, CXn2, CXn3 and CXi determine the saturation limit, half maximum, slope and initial value of the CX 1 curve, respectively; and CXn4 is the parameter for the time constant of CX dynamics.

Parameter setting
In this modeling framework, the first module (for calcium dynamics) and the third module (for force production) were constrained by the relevant experimental data. Only the second module (representing the calcium-force relationship) was adjusted for the model to replicate the force data measured under physiological conditions. First, the parameters of module 3 were determined from the data for length-and velocity-tension properties under full excitation. The K SE for the stiffness of serial elastic elements was estimated by dividing the difference in force (i.e., 30.6 N for CAT14 and 9.5 N for CAT12) by the variation in muscle-tendon length (i.e., 1.9 mm for CAT14 and 1.55 mm for CAT12) and then normalized with the peak isometric force (P 0.5 ) at X m,0.5 (i.e., 100.3 N for CAT14 and 17.5 N for CAT12).  16), respectively, in [21]). Second, all parameter values of module 1 except τ 1 and τ 2 for the rate of calcium release from the sarcoplasmic reticulum were adopted from the previous study in [21]. τ 1 and τ 2 of module 1 and C1i, C2i, C3-C5 and α of module 2 were simultaneously optimized to fit the shape of the twitch and tetanic responses during isometric contraction at X m,0.5 for stimulation frequencies of 10 and 40 Hz for CAT14 and 20 and 40 Hz for CAT12, without considering sag behavior by setting C1n1 and C2n1 to 0. The force data obtained at 20 Hz for CAT14 and 10 Hz for CAT12 were not involved because the muscle showed the strongest sag behavior under this stimulation frequency. Third, C1n1-C1n4 in module 2 were first optimized to replicate the sag behavior observed during the tetanic response at 20 Hz for CAT14 and 10 Hz for CAT12, and then C2n1-C2n4 were further optimized for the tetanic response at 40 Hz for both CAT14 and CAT12 under isometric conditions at X m,0.5 . The optimization of parameter values was performed using the optimization tool (i.e., Multiple Run Fitter based on the PRAXIS algorithm) built in NEURON software [39]. Finally, φ 1 -φ 4 of module 1 for the length-dependent changes in the calcium-force relationship were determined to match the amplitude of the twitch at the physiologically minimal, intermediate and maximal muscle lengths using the polyfit function based on the least squares method built in MATLAB software. Notably, the model parameters (i.e., g1-g3 and a0-d0) in Module 3 were forced to reproduce the length-tension (L-T) and velocity-tension (V-T) relationship measured under full excitation with 100 Hz frequency stimulation. From the perspective of sliding filament theory [40,41], the values of those parameters reflected the number (N) of force-generating crossbridges formed from spatiotemporal actin-myosin interactions underlying the L-T and V-T relationships. However, for partial excitation with low-frequency stimulation, changes in the fraction of N were reflected by the A-CaT relationship of Module 2. In this study, thus, the number of crossbridges in the force-generating state at a given time corresponds to the product of the maximal tension produced in a given length and velocity (T(L, V)) and A under physiological input conditions.
All parameter values used for the model of the medial gastrocnemius muscle are presented in Table 1.

Comparison of experimental and simulation data
For the purpose of comparison between different cat MG muscles, all active forces measured experimentally were normalized with P 0.5 [32]. The error between the experimental and simulation data was quantitatively compared across different levels of current stimulation by calculating the root mean square error normalized with maximal force generation using the where the subscripts m and s and N indicate measurement, simulation, and the number of data points (indexed with i) considered for evaluation, respectively.

Numerical approaches
The model of fast muscle was implemented using the NEURON model description language (NMODL) in the NEURON software environment (version 6.1). The numerical integration of model equations was performed using the integration method (cnexp) built in NEURON with a fixed time step of 0.025 ms. Under this simulation condition, the stability and accuracy of simulations were confirmed while varying parameter values over a wide range. The initial values for Ca SR and CS in the SR and B and T in the SP were set to 2.5 mM, 30 mM, 0.43 mM, and 70 μM, respectively, based on a previous experimental study [33]. The codes for simulations of the muscle-tendon model developed in this study are presented as supplemental information (S1 Data) and also available in the public repository of ModelDB (http://modeldb.yale. edu/267738).

Results
We constructed a computational muscle-tendon model comprising three submodules to demonstrate whether and how dynamic alterations in the calcium-force relationship of the myofilaments could explain force variations experimentally observed in fast skeletal muscles under physiological conditions (Fig 1). The calcium-force relationship experimentally characterized in skinned muscle fibers was represented as the steady-state relationship of the muscle activation level (i.e.,Ã 1 ) to the constant concentration of calcium-binding troponin (i.e., CaT/T 0 ) in the modular model for skeletal muscle (module 2 in Fig 1A). The CaT/T 0 for half-maximal activation (i.e., C1 comparable to pCa in the Hill function) and the slope at the CaT/T 0 for half-maximal activation (i.e., 1/C2 comparable to Hill coefficient, n H ) in theÃ 1 -CaT/T 0 relationship were dynamically varied to infer the dynamic alterations in calcium sensitivity and cooperativity during force production under physiological conditions. The C1 and C2 variations were modulated as a function of CaT/T 0 considering the saturation (i.e., C1n1 and C2n1), threshold (i.e., C1n2 and C2n2), and rate (i.e., τ C1 and τ C2 ) in the modular muscle-tendon model (see Fig 1B for graphical illustration and the Methods for mathematical equations).

Absence of sag behavior under a static calcium-force relationship
First, it was evaluated whether the model could reproduce the muscle force obtained during isometric contractions at the intermediate length (i.e., X m,0.5 = 5 mm) between the physiological maximum (i.e., X m,1 = 10 mm) and minimum (i.e., X m,0 = 0 mm) over a full physiological range of stimulation frequencies (i.e., 1~40 Hz) with no change in the calcium-force relationship (Fig 2). This condition was implemented by setting the values of C1n1 and C2n1 to zero. The parameters (i.e., τ 1 and τ 2 of module 1 and C1i, C2i, C3, C4, C5 and α of module 2) shaping the static calcium-force relationship were optimized to match the experimental data showing no sag phenomenon obtained at 1, 10, and 40 Hz from the medial gastrocnemius of a cat. The model with the static calcium-force relationship could predict the time course of the force profile measured at 1, 10 and 40 as well as 100 Hz. However, the error between the experimental and simulated forces was maximized at a stimulation frequency of 20 Hz, showing the absence of the sag phenomenon. This result indicates the insufficiency of the previous model, with no variation in the calcium-force relation, for predicting force production by fast muscle during unfused isometric contractions.

Depression of calcium sensitivity under low excitation
Then, we assessed whether the error between the experiment and simulation at 20 Hz could be compensated by dynamically varying the calcium-force relationship. The parameters (i.e., C1n1-

PLOS COMPUTATIONAL BIOLOGY
Dynamic calcium-force relationship in fast muscle C1n4) underlying the dynamic variation of C1 were optimized for the force data at the stimulation frequency of 20 Hz. The simulation error at 20 Hz could be dramatically reduced by progressively increasing C1 (Fig 3). The variation in 1/C2 alone through parameter optimization (i.e., C2n1-C2n4) for the force data at 20 Hz was not sufficient to reproduce the progressive force decline observed during unfused tetanic contractions at 20 Hz (S1 Fig). The decrease in calcium sensitivity, rather than cooperativity, could more accurately predict force production under low muscle excitation levels (i.e., 20 Hz). The dynamic increase in C1 did not affect the force production at stimulation frequencies less than 20 Hz. However, the simulation error was significantly increased under higher levels of stimulation frequency (i.e., > 20 Hz). The dynamic increase in C1 alone was not sufficient to reproduce the time course of force production at 40 Hz.

Potentiation of calcium cooperativity under high excitation
Finally, we assessed whether the increased error between the experiment and simulation with only variation in C1 for stimulation frequencies > 20 Hz could be compensated by dynamically varying the calcium-force relationship. The parameters (i.e., C2n1-C2n4) underlying the dynamic variation of 1/C2 were further optimized for the force data at the stimulation frequency of 40 Hz. The simulation error with only variation in C1 for stimulation frequencies > 20 Hz, including 100 Hz, could be substantially reduced by dynamically increasing the 1/C2 value (Fig 4). The variation in C1 alone through parameter optimization (i.e., C1n1-C1n4) for the force data at 40 Hz was not sufficient to reproduce the slow rate of force development following the fast rate of initial force development at 40 Hz (S2 Fig). The gradual enhancement of calcium cooperativity, rather than sensitivity, could lead to better prediction of the slow force development observed at 40 Hz. In addition, the forces at lower stimulation frequencies < 40 Hz could be simulated without the slope variation of the calcium-force relationship by making its threshold a function of the calcium-bound troponin (i.e., CaT) concentration. All these results indicate that calcium sensitivity and cooperativity operationally vary depending upon the state of muscle excitation under physiological conditions.

Dynamics of the calcium-force relationship at the intermediate length
To visualize how the calcium-force relationship varies depending on the level of excitation, the steady-state calcium-activation (i.e., CaT/T 0 -Ã 1 ) curve was reconstructed at the initial and maximally varied states during isometric contractions at the intermediate muscle-tendon length over the full physiological range of stimulation frequencies (Fig 5A). The shift of the midpoint of the calcium-activation relation to the right was facilitated as the excitation level increased from 20 Hz. However, the slope of the calcium-activation relationship was steepened at higher levels of excitation frequencies from 40 Hz. These findings obtained under steadystate conditions were also confirmed in the transient state for the calcium binding troponinmuscle activation (i.e., CaT/T 0 -A) ( Fig 5B) and intracellular calcium-muscle force (i.e., Ca-F) (Fig 5C) relationship. The steady-state calcium-activation relationships could be reflected in the transient calcium-force relationships measured during the relaxation of force production (see the insets in Fig 5C). These results imply that the increase in the threshold for crossbridge formation can be counteracted by an enhancement in the cooperativity for crossbridge formation as the level of muscle excitation increases.

Role of calcium cooperativity in shaping sag behavior
Compared to the experimental data obtained from the cat medial gastrocnemius, we have shown the capability of the muscle-tendon model to produce a typical form of sag behavior

PLOS COMPUTATIONAL BIOLOGY
Dynamic calcium-force relationship in fast muscle showing progressive force decline after the initial peak during the unfused isometric contraction. It was further assessed whether the model could represent more complex forms of sag behavior that have been reported in previous experimental studies [42,43]. Fig 6 shows the four representative types of sag behavior that are identifiable by the model during unfused

PLOS COMPUTATIONAL BIOLOGY
Dynamic calcium-force relationship in fast muscle isometric contractions at the intermediate length for stimulation frequencies ranging from 20 to 30 Hz. The simple form (referred to as Type I) of sag behavior with progressive force decline after the initial peak could be produced by changing C1 (Fig 6A). However, complex forms of sag behavior could be induced by varying the C2 parameters at the same time. With a lower threshold (less C2n2) and slower speed (larger τ C2 ) (i.e., 100 ms) in the slope variation of the calcium-force relationship, the simple sag behavior transitioned to a complex form (referred to as Type II) of sag behavior, showing multiple force peaks following the initial force peak at 20 Hz (Fig 6B). With the same C2n2 but much slower τ C2 (i.e., 500 ms) in comparison to the Type II sag, another complex form (referred to as Type III) of sag behavior showing double force declines following the initial force peak was induced at 20 Hz (Fig 6C). With the same values of C2n2 and τ C2 used for the Type III sag, the last complex form (referred to as Type IV) of force decline followed by force amplification after the initial force peak was determined to be under 30 Hz (Fig 6D). Type I, II and IV sags have been experimentally identified. In contrast, Type III sag theoretically emerged through the analysis of the muscle-tendon model in this study. We further investigated how C2n2, τ C2 and excitation level might interact to shape the sag behavior during submaximal contractions at the intermediate muscle-tendon length. Fig 7 shows the map of sag types over the parameter space for C2n2 and τ C2 under two different levels (i.e., 20 and 30 Hz) of muscle excitation. Type I sag tended to transition to Type II sag as C2n2 and τ C2 decreased, while Type I sag tended to transition to Type III sag as C2n2 decreased and τ C2 increased under a stimulation frequency of 20 Hz. Interestingly, the interaction of C2n2 and τ C2 on sag type differed under the stimulation frequency of 30 Hz. The sag phenomenon tended to disappear as C2n2 and τ C2 decreased, while Type I sag tended to transition to Type IV as C2n2 decreased and τ C2 increased. These results indicate that the type of sag phenomenon might be determined predominantly by a dynamic change in the slope of the calcium-force relationship depending on the level of excitation.

Calcium-force relationship for the L-T and V-T properties
The dynamic variations in the calcium-force relationship were further estimated for the force production under full excitation (i.e., 100 Hz) while varying the muscle-tendon length at a constant velocity. The muscle-tendon model with variations in the calcium-force relationship could accurately demonstrate the characteristics of force-length and force-velocity properties, as seen in the experiment under full excitation (Fig 8). This result was attributed to the fact that the depression effect of the rightward drift (i.e., C1) of the calcium-force relationship (S3 Fig) was sufficiently compensated by the potentiation effect of the upward drift (i.e., C2) of the slope in the calcium-force relationship under full excitation. As a result, the dynamic variation in the calcium-force relationship did not tend to affect the force production because the activation level was almost maximal under full excitation. This result suggests the importance of the net effect of both changes in the midpoint and the slope of the calcium-force relationship in predicting the length-and velocity-tension properties of fast skeletal muscle under physiological conditions.

Model predictions for a different cat
We tested whether the results hold for cat MG muscles showing different sag behavior during unfused isometric contractions. The same modeling approach used for CAT14 was applied to the force data independently measured from the different cat MG muscles (i.e., CAT12) with the sag behavior occurring at the lower level of stimulation frequency (i.e., 10 Hz) than CAT14 (see Table 1 for the parameter values of model CAT12). As expected from the model of CAT14, the muscle-tendon model for the new force data could not exhibit the sag behavior without variations in the calcium-force relationship during isometric contractions at the intermediate muscle-tendon length (i.e., X m,0.5 ). The parameters for C1 dynamics (i.e., C1n1-C1n4) and C2 dynamics (i.e., C2n1-C2n4) were optimized to best represent the new force data at stimulation frequencies of 10 and 40 Hz, respectively. The midpoint (i.e., C1) of the calciumforce relationship was required to drift rightward, whereas the slope (i.e., 1/C2) of the calciumforce relationship needed to drift upward to reproduce the isometric forces over a full range of stimulation frequencies (i.e., 1-100 Hz) at X m,0.5 , including sag behavior (S4 Fig). The pattern of dynamic alterations in the calcium-force relationship was similar to that inferred from the dataset for CAT14 (S5 Fig). All types of sag behavior identified by model CAT14 were also found by varying the parameters (i.e., C2n2 and τ C2 ) for the slope variation of the calciumforce relationship at different excitation levels (i.e., 20 and 30 Hz) in model CAT12 (S6 Fig). In addition, the model CAT12 reflecting the length dependent twitch response could also replicate the force-length and force-velocity properties under full excitation, showing the importance of the net effect of C1 and C2 variation in the calcium-force relationship (S7 Fig). To this end, the characteristics of dynamic alterations in the calcium-force relationship seem to be similar across fast-type muscles with sag behavior in adult cats.

PLOS COMPUTATIONAL BIOLOGY
Dynamic calcium-force relationship in fast muscle

Dynamic calcium-force relationship under physiological conditions
The dynamics of the calcium-force relationship during force production differed depending upon the muscle type (slow-twitch versus fast-twitch). In contrast to the case for slow muscles such as the soleus [21], our modular model of the cat medial gastrocnemius muscle predicted that the calcium sensitivity (assessed by C1) and cooperativity (assessed by 1/C2) should dynamically change during force generation over a physiological range of stimulation frequencies (20~40 Hz) and muscle-tendon lengths (0~10 mm) in fast skeletal muscle.
The static calcium-force relationship has been analyzed to infer the calcium dependence of force-producing crossbridge formation under in vitro steady-state conditions [44]. Our study may provide insights into the dynamic aspects regarding the calcium sensitivity and cooperativity of crossbridge formation under physiological conditions. That is, the calcium sensitivity and cooperativity of crossbridge formation might tend to be constant under low levels of muscle excitation (i.e., < 10 Hz) (Fig 2). However, the calcium sensitivity of crossbridge formation might be progressively lowered, resulting in sag behavior under intermediate levels of muscle excitation (i.e., approximately 20 Hz) (Fig 3). The cooperativity of crossbridge formation might be progressively increased, resulting in force enhancement under high levels of muscle excitation (i.e., > 40 Hz) during isometric and isokinetic contractions (Figs 4 and 8). Furthermore, the cooperativity of crossbridge formation might play a crucial role in determining the form of sag behavior (Figs 6 and 7).
All these findings suggest that the calcium sensitivity and cooperativity of crossbridge formation might operationally vary depending on both neural excitation and muscle length in fast muscle under physiological conditions.

Factors regulating the dynamics of the calcium-force relationship
Contractile force is thought to be determined by the number and force of strongly bound crossbridges and the rate of crossbridge cycling. These crossbridge properties are likely to be regulated dynamically by multiple factors. The rightward drift of the calcium-force relationship might result from a gradual reduction in the number of force-generating crossbridges and the force per force-generating crossbridge under low pH or high inorganic phosphate concentration [45] or in the duration of strong actomyosin binding during crossbridge cycling [19]. The upward drift of the slope of the calcium-force relationship might occur due to a progressive increase in the transition rate from non-force-generating to force-generating crossbridge formation through myosin light chain phosphorylation [46] or in the binding affinity of myosin heads to actin filaments due to mechanosensitive binding proteins [47]. Slow-and fasttwitch muscle fibers have been shown to correspond to Type I and Type II chemically identified by myosin ATPase or MHC I and MHC II by myosin heavy chain isoform [48]. The typespecific compositions of myosin proteins affecting crossbridge cycling have been suggested to underlie slow but sustained contractions of slow muscle fibers and fast but fatiguing contractions of fast muscle fibers [49,50]. These findings from isolated muscle fibers might also explain the differential dynamics of the calcium-force relationship between the fast (Fig 8) and slow (Fig 3 in [21]) muscle models for the length-and velocity-tension relationships under full excitation. The modular modeling framework developed for this study may provide a testbed for further investigation of the relative or combinational contributions of these mechanisms to the calcium-force relationship under physiological conditions.

Comparison with previous studies
Previous experimental studies have shown a rightward midpoint shift and a decreased slope steepness in the calcium-force relationship as the inorganic phosphate (P i ) concentration increases in isolated muscle fibers [18]. As shown in the previous study, the midpoint (i.e., C1) of the calcium-force relationship was required to shift to the right as the excitation level increased in the present study. In contrast to the previous study, however, the slope (i.e., 1/C2) of the calcium-force relationship needed to be steepened as the excitation level increased under physiological conditions (Fig 4). This result suggests that the cooperativity of crossbridge formation might be enhanced through other factors, such as strain-induced potentiation [51,52], counteracting the inhibitory effects of P i accumulation at high excitation levels in fast muscle under physiological conditions. Sag behavior has been investigated in fast muscles in cats [22], rats [42], frogs [53] and humans [54]. A variety of sag forms have been identified from rat triceps surae motor units, including one simple form and two complex forms of the sag phenomenon [42,43]. In this study, our model could reproduce all three forms of the sag phenomenon, suggesting that the simple form (i.e., Type I) tended to be determined by a rightward shift (i.e., C1) of the calcium-force relationship, whereas the complex forms (i.e., Type II & Iv) tended to be shaped by the slope variation (i.e., 1/C2) in the calcium-force relationship (Fig 6). The muscle-tendon model could further identify different types (e.g., Type III) of sag behavior while varying the excitation level and the parameters (i.e., C2n2 and τ C2 ) for the threshold and speed in the slope variation of the calcium-force relationship (Fig 7). This result suggests that the interplay between the excitation level and the slope variation in the calcium-force relationship plays a critical role in determining the form of the sag phenomenon.
In intact animal preparations, the sag magnitude has been shown to be greater at shortened lengths than at lengthened lengths [55]. In isolated muscle fiber preparations, the slope of the steady-state calcium-force relationship has been reported to decrease as the muscle length is stretched [7]. Thus, we tested whether the length dependence of sag magnitude could be explained by the slope variation of the calcium-force relationship during isometric contractions at 20 Hz for CAT14 and 10 Hz for CAT12. The model without slope variation in the calcium-force relationship exhibited the increasing sag amplitude with increasing muscle length during isomeric contractions at 20 Hz (Figs 9A1 and 9B1 for CAT14 and S8A1 and S8B1 for CAT12). When the slope of the calcium-force relationship was varied by lowering C2n2 and decreasing C2n1 in the shortened state and increasing C2n1 in the lengthened state from the default value for the intermediate state, as reported from skinned fiber experiments [7], the model muscle displayed the opposite tendency, showing a decreasing sag amplitude with increasing muscle length (Figs 9A2 and 9B2 for CAT14 and S8A2 and S8B2 for CAT12). The length-dependent sag magnitude could be explained by the progressively increased calcium cooperativity at shortened muscle lengths but progressively lowered calcium cooperativity at lengthened muscle lengths under intermediate levels of muscle excitation. This result highlighted the interplay between the slope of the calcium-force relationship and muscle length in shaping the sag behavior under physiological conditions. The computational model for fast skeletal muscle has been developed phenomenologically, focusing on capturing force-length and force-velocity properties across different levels of muscle excitation [55] or superposition of twitch responses under isometric conditions [56,57]. More recently, it has been proposed that the cooperativity of crossbridge formation should be incorporated into the phenomenological model to reflect the physiological conditions of muscle activation dynamics [58]. In the present study, we have further demonstrated that the dynamic variation in the calcium-force relationship must be considered to reproduce the complex forms of force production experimentally shown in fast muscle studies under submaximal rather than maximal contractions (see Fig 8 for the force-length and force-velocity properties under full excitation). Therefore, the muscle-tendon model developed for this study may provide a physiologically plausible framework for force production over a full physiological range of stimulation frequencies and muscle lengths under isometric and isokinetic conditions.

Influence of Modules 1 and 3 on sag behavior
In Module 1, we employed a canonical model to transform electrical stimulation into the concentration of calcium binding troponin. This model assumed that calcium ions rapidly occupy the regulatory sites of troponin, and the total troponin concentration is similar in both slowand fast-twitch muscle fibers [59]. However, type-specific properties, such as the dominant existence of parvalbumin and two regulatory sites of troponin in fast muscle fibers [60], might play a role in determining the sag phenomenon.
We first evaluated the influence of parvalbumin with no variation in C1 and C2 (C1n1 = C2n1 = 0). The simulations in Fig 2 were repeated with the same kinetic scheme (K3 = 4.17*10 7 M -1 *ms -1 and K4 = 0.5 ms -1 ) and total concentration (B 0 = 7.5*10 −4 M) as reported in a previous study [60]. The results showed a substantial reduction (~40% on average) in predicted peak forces, particularly for low stimulation frequencies (< 40 Hz), without inducing sag behavior (S9 Fig). Then, we assessed the influence of two regulatory sites of troponin. The reaction equation of Ca + T = CaT was changed to 2Ca + T = Ca 2 T in Module 1 to explicitly express the interaction between free calcium ions and two troponin regulatory sites ( Fig 1A). Accordingly, the system Eqs (3) and (5) were updated as follows, The total troponin concentration (T 0 = 1.2*10 −4 M) was set to be the same as that used in a previous study [60]. The twitch response was matched to the force data measured at X m,0.5 by increasing the calcium release ([Ca SR ] init = 0.04 M, R max = 139 ms -1 and τ 2 = 14.5 ms) from the SR to compensate for the increase in parvalbumin and troponin site concentration. The reaction constants of K5i and K6i were reset to 1.37*10 6 M -2 ms -1 and 0.56 ms -1 so that the calcium concentration for half activation (Ca 50 = sqrt(K6i/K5i) = 6.4*10 −4 M) became 70% larger and the calcium-force relationship slope steeper than the canonical case (Ca 50 = K6i/ K5i = 3.75*10 −4 M). The degree of Ca 50 increase was determined as the average increase in Ca 50 in fast muscle fibers compared to slow muscle fibers across rat (100% increase in [61]), monkey (35% increase in [62]), and human (80% increase in [17]) skeletal muscle. As a result, the new model with fast-type-specific properties could not reproduce the sag behavior without variations in C1 and C2 and produced larger forces at high stimulation frequencies (> 20 Hz) due to the steeper calcium-force relationship (S10 Fig). However, with C1i (0.162), C2i (0.134), C2n1 (-0.0095), and C2n4 (200 ms) in Module 2 modestly adjusted, the new fast muscle model could reproduce all isometric and isokinetic force data showing dynamic variations in C1 and C2 similar to those for the canonical case (S11 and S12 Figs).
In addition, calcium binding (K5) and unbinding (K6) to troponin were represented as a function of muscle-tendon length (X m ) and muscle activation level (Ã) in Module 1. The isometric forces in Fig 2 were simulated at shortened (X m,0 ) and lengthened (X m,1 ) lengths with and without the length dependence of K5 (K5i for all X m ) under no variation in C1 and C2 (S13 Fig). We also simulated the isometric forces at the intermediate length (X m,0.5 ) with and without the activation dependence of K6 (K6i for allÃ) without C1 and C2 variation (S14 Fig). The length dependence of K5 had a modest impact on force predictions only for low stimulation frequencies (< 40 Hz) at the shortened state without inducing the sag behavior. However, the lack of activation dependence in K6 caused a considerable reduction (~85 ± 0.08%, mean ± SD) in predicted peak forces for all stimulation frequencies without producing the sag behavior.
Furthermore, the present study did not explicitly consider titin-myofilament interactions for predicting active force production during isometric contractions. Instead, extracellular (tendon and aponeurosis) and intracellular (myofilaments and titin proteins) elastic properties were collectively reflected in the serial elastic stiffness (K se ) in Module 3. In the present study, K se was constant but may vary partially due to changes in titin stiffness upon muscle activation [63]. All simulations in Figs 4 and 8 were repeated over a range of K se values (0.13~0.23 mm -1 ) adopted from a previous study measuring the short-range stiffness as a function of force for cat medial gastrocnemius muscles [64]. The variation in K se over the physiological range had little impact on the predictions of forces during isometric and isokinetic contractions. Notably, titin-myofilament interactions have not affected the length-dependent calcium activation and force-velocity properties [65]. Titin has also been suggested as an active molecular spring underlying the activity-dependent response of the muscle, such as residual force enhancement or suppression after active stretching or shortening [66]. This issue is beyond the scope of the present study focusing on active force production under isometric and isokinetic conditions.
Regarding Module 3, the parallel passive elements were not included in this study. The reason was that the passive force measured without excitation was very low and constant compared to the active force over our investigation's range of muscle length (S15A Fig). In addition, the sag profile, including both active and passive force, was shifted downward when compared with the active force at the same muscle-tendon length (S15B Fig). This result indicates that the dampening and velocity dependency of parallel elements may not significantly impact sag behavior.
All these results reinforce the hypothesis that both calcium sensitivity and cooperativity of crossbridge formation dynamically change to produce the sag behavior observed in fast skeletal muscles under physiological conditions.

Current limitations
The model including the dynamic variation in the calcium-force relationship could reproduce the force production of a cat medial gastrocnemius muscle over a full physiological range of stimulation frequencies and muscle lengths under isometric and isokinetic conditions. However, this model should be further improved to reflect multiple factors determining the physiological behaviors of fast muscles. First, in the present study, the mechanisms controlling intracellular calcium concentration (Module 1) were assumed not to change during submaximal muscle contractions over short periods of time (i.e., < 1 second). However, this assumption may not hold for maximal muscle contractions over longer periods of time (i.e., > several seconds) since under this condition, the intracellular calcium concentration has been reported to increase and then decrease during progressive force decline or fatigue [24]. Second, the calcium-force relationship was dynamically varied by directly manipulating the parameters of the curve (Module 2), which simulated a progressive increase in the threshold and the cooperativity for actomyosin crossbridge formation. However, the modular modeling approach used for this study could be extended to systematically evaluate potential subcellular mechanisms that might underlie the dynamic variation in the calcium-force relationship during force production in fast muscles [67]. Third, the calcium dynamics in the sarcoplasm were modeled based on the experimental data obtained from rats in the present model. Further experimental investigations are needed to identify the specific differences in sarcoplasmic calcium dynamics between different species and muscle types [68]. Fourth, the modeling approach used for this study represents the average behavior of whole muscle comprising different muscle fibers.
Although cat gastrocnemius muscles are known to be composed predominantly of fast-twitch muscle fibers (i.e., over 75%) [31], this limitation might increase the error between the simulation and experimental data. Future studies should be performed to determine how the heterogeneous organization of whole muscle influences the force output [20]. Fifth, the current model was validated only for constant stimulation frequencies. Neural excitation has been shown to be random rather than constant [69]. Thus, the current model might not reflect the random effect of neural excitation on force production, such as the catch-like property [70]. Finally, the present model is limited to force production under isometric and isokinetic muscle conditions. Future experiments are needed to measure the data on force production during dynamic variations in muscle length, such as locomotor-like movements [29].

Conclusions
Unlike in slow-twitch muscles, our model suggests that the calcium-force relationship dynamically changes during the force responses of fast-twitch muscles over a full physiological range of stimulation frequencies and muscle lengths. The dynamic variation in the calcium-force relationship operationally varied depending upon the excitation level and muscle length during unfused isometric contractions, whereas it did not tend to affect the length-and velocitytension properties under full excitation. The slope variation in the calcium-force relationship was crucial to determine the types and length dependence of sag phenomena. These results provide insights into the dynamic properties of crossbridge formation during force production and a computational framework for physiologically realistic modeling of skeletal muscles.